#include "hamiltonian.h"

using arma::cx_vec;

// This operator (use by simply appending () to a Hamiltonian
// type variable) calculates the transfer integral matrix H 
// as a sparse complex matrix (sp_cx_mat). The operator returns
// the object it self for easy access to result.
arma_noinline Hamiltonian &Hamiltonian::operator()(int id)
{
	// First we get the unit from the id.
	Unit* unitref = system->getUnit(id);
	if(!unitref) // If the unit is invalid then we do nothing.
		return *this;
	
	// Now we make a vector for the values of the hamilton matrix.
	const size_t elementsize = unitref->getBonds().size();
	cx_vec values = arma::zeros<cx_vec>(elementsize); // All values are initialized to zero.

	// We iterate (parallel) across the elements, and add values for each bond in the given element.
	_Cilk_for(size_t e = 0; e < elementsize; e++)
	{
		// For each element there is a set of bonds. The final element is a sum of these bonds.
		size_t bondsize = unitref->getBonds()[e].size();
		for(size_t b = 0; b < bondsize; b++)
			values(e) += getHamiltonElement(unitref->getBonds()[e][b]);
	}

	// We set the result using precalculated locations (i,j corresponding to the values).
	H = sp_cx_mat(unitref->getBondsLocations(),values,false);

	return *this;
}

Hamiltonian::Hamiltonian() : 
	// A pointer to the used system. The system contains variables 
	// needed to construct elements (coupling strengths etc.).
	system(&System::Std),
	// A double containing the current bias (for bi-layered 
	// systems).
	bias(0.),
	// A double vector containing the current k-point.
	k(arma::zeros<vec>(3)) 
{ }
	
complex Hamiltonian::getHamiltonElement(const Bond &bond)
{
	// First get the phase phi from the k-point,
	double phi = arma::dot(k,bond.xyz());
	// then translate to the complex constant c.
	complex c(cos(phi),sin(phi));

	if(bond.isBetweenSameAtom()) // The bond origin and destination atoms are identical.
	{
		// The energy is shifte according to bias and dimer shift (Delta).
		// Note that odd index layers (top of bi-layer) have positive onsite shift,
		// while even index layers (bottom of bi-layer) have negative onsite shift.
		// h() is the height index.
		if(bond.getOrigin()->isDimer())
			return system->epsilon + system->Delta + ((bond.getOrigin()->h() % 2) ? bias/2. : -bias/2. );
		else
			return system->epsilon + ((bond.getOrigin()->h() % 2) ? bias/2. : -bias/2. );
	}
	else if(bond.isBetweenInplaneNN())
		return system->gamma0*c;
	else if(bond.isBetweenDimerNN())
		return system->gamma1*c;
	else if(bond.isBetweenNonDimerNN())
		return system->gamma3*c;
	else if(bond.isBetweenDimerNonDimerNN())
		return system->gamma4*c;

	return 0.;
}